# Loughney et al., "Tectonic influence on Cenozoic mammal richness and sedimentation history of the Basin and Range, western North America"
# Science Advances 
# R functions file

# 1. richnessBoot() randomly samples the age of fossil localities and tallies the number of species occurrences in each time bin
# arguments 
# - occurrences: .csv file of species occurrences; 
# - numBins: number of time bins; 
# - binSize: time-bin duration (default is 0.5 Myr)
# returns
# - divTrial: matrix of number of species in each time bin
richnessBoot <- function(occurrences, numBins, binSize = 0.5) {
	locAge <- runif(nrow(occurrences), occurrences$MinAge, occurrences$MaxAge)
	occurrences$trialAge <- locAge
	bin <- ceiling(occurrences$trialAge / binSize)
	youngestAge <- bin * binSize - binSize
	oldestAge <- bin * binSize
	assBins <- ifelse(occurrences$trialAge <= oldestAge, floor(occurrences$trialAge / binSize) + 1, floor(occurrences$trialAge / binSize))
	occurrences$bin <- assBins
	divTrial <- rep(0, numBins)	
	for(i in 1:numBins) {
		divTrial[i] <- length(unique(occurrences$SpName[which(occurrences$bin == i )]))
		}
	return(divTrial)	
	}

# 2. localityBoot() randomly samples age of fossil localities and tallies the number of localities in each time bin
# arguments 
# - occurrences: .csv file of species occurrences; 
# - numBins: number of time bins; 
# - binSize: time-bin duration (default is 0.5 Myr)
# returns
# - locTrial: matrix of number of localities in each time bin
localityBoot <- function(occurrences, numBins, binSize = 0.5) {
	locAge <- runif(nrow(occurrences), occurrences$MinAge, occurrences$MaxAge)
	occurrences$trialAge <- locAge
	bin <- ceiling(occurrences$trialAge / binSize)
	youngestAge <- bin * binSize - binSize
	oldestAge <- bin * binSize
	assBins <- ifelse(occurrences$trialAge <= oldestAge, floor(occurrences$trialAge / binSize) + 1, floor(occurrences$trialAge / binSize))
	occurrences$bin <- assBins
	locTrial <- rep(0, numBins)	
	for(i in 1:numBins) {
		locTrial[i] <- length(unique(occurrences$Locality[which(occurrences$bin == i )]))
		}
	return(locTrial)	
	}

# 3. ciQuantiles() calculates quantiles on number of species richness or localities in each time bin
# arguments 
# - diversity: bootstrapped species richness or localities; 
# - confidence: the confidence level; 
# returns
# - divCurves: data frame of 2.5% richness, 50% richness, 97.5% richness, age
ciQuantiles <- function(diversity, confidence) {
	CIup <- confidence + (1 - confidence)/2
	CIlo <- (1 - confidence)/2
	quantiles <- c(CIlo, 0.5, CIup)
	divQuantiles <- t(apply(diversity, MARGIN = 2, quantile, probs = quantiles))
	divCurves <- as.data.frame(divQuantiles)
	topAges <- seq(0, 35.5, 0.5)
	divCurves$age <- topAges
	return(divCurves)
	}
